Dynamics of quantum dissipation systems interacting with Fermion and Boson grand 
canonical bath ensembles: Hierarchical equations of motion approach 
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A hierarchical equations of motion formalism for a quantum dissipation system in a grand canon- 
ical bath ensemble surrounding is constructed, on the basis of the calculus-on-path-integral algo- 
rithm, together with the parametrization of arbitrary non-Markovin bath that satisfies fluctuation- 
dissipation theorem. The influence functionals for both the Fermion or Boson bath interaction are 
found to be of the same path-integral expression as the canonical bath, assuming they all satisfy the 
Gaussian statistics. However, the equation of motion formalism are different, due to the fluctuation- 
dissipation theories that are distinct and used explicitly. The implications of the present work to 
quantum transport through molecular wires and electron transfer in complex molecular systems are 
discussed. 
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I. INTRODUCTION 

Quantum dissipation theory (QDT) is concerned with 
the fundamental formulations for the dynamics of a quan- 
tum system of primary interest embedded in a quan- 
tum bath surrounding environment. The key quantity 
in QDT is the so-called reduced system density op- 
erator, p{t) = tr B /5x(i) ; ke., the partial trace of the 
total density operator over the bath space of practi- 
cally infinite degrees of freedom. Due to its fundamen- 
tal importance and intrinsic complexity, the develop- 
ment of QDT has involved scientists from many fields 
of research J ii2 1 3 1 4 1 5 1 6 i 7 1 8 1 9 i io ^ j-^g rem ained as a chal- 
lenging topic since the middle of the last century. 

The exact p(t) can be expressed in terms of the 
Feynman- Vernon influence functional, assuming the in- 
teraction bath satisfies the Gaussian statistics ^ Its nu- 
merical implementation has been carried out in a few 
simple systems with the forward-backward iterative path 
integral propagation methodi 12 i 13 Unraveling the path 
integral formalism into stochastic Schrodinger or stochas- 
tic Liouville equations has also been proposed in various 
forms in the past 10 years, 14 i 15 i 16 i 17 i 18 i 19 i 20 i 21 i 22 i 23 ' 24 i 25 
The main obstacle to both the path integral and stochas- 
tic differential equation QDTs is their formidable numer- 
ical implementation in multilevel systems. Applications 
based on the path integral or stochastic QDT formalism 
are also often too cumbersome to be practical. 

For the numerical efficiency and practical use in 
general, one will be interested in a differential equa- 
tions of motion (EOM) formalism. In principle, a for- 
mally exact (but unclosed) EOM formalism of QDT 
can be constructed via the Nakajima-Zwanzig-Mori pro- 
jection operator technique ; 26 i 27 ' 28 i 29 To complete the 
formalism, however, one usually has to invoke cer- 
tain approximations. In particular, various commonly 
used forms of Bloch-Redfield theory and Fokker-Planck 
equation o 30 ! 31 ^^^ 4 .! 3 ^^ involve the Markovian limit 



and the weak system-bath interaction. The validi- 
ties of these approximations are increasingly challenged, 
especially due to the emerging fields of nanoscience 
such as quantum transport and quantum information 
procc::c: . 37,38 - 39 - 40 - 41 - 42 - 43,44,45,46 - 47,48,49 

We have recently constructed a formally exact 
EOM formalism for arbitrary non-Markovian dissipa- 
tion systems interacting with Gaussian canonical bath 
ensemblei 50 ! 51 This theory generalizes the Tanimura and 
coworkers' hierarchical formalism i 52 i 53 ' 54 i 55 The present 
paper is to extend our previous results to include both 
the Fermion and Boson grand canonical bath ensemble 
cases. The desired hierarchical EOM will be constructed, 
via the auxiliary influence-generating functionals (IGF) 
approach, together with the calculus-on-path-integrals 
(COPI) algorithm.^ 1 - 

The remainder of this paper is organized as fol- 
lows. In Sec. HIl after the general description of a 
system-bath composite Hamiltonians, we discuss the 
fluctuation-dissipation theorem and other useful quan- 
tum mechanics relations for the grand canonical ensem- 
bles of Fermion/Boson particles. In Sec. IIIIl we re- 
visit the influence functionals path-integral formalism of 
QDT, assuming that the interaction grand canonical bath 
satisfies Gaussian statistics. The derivations in relation 
to these two sections are detailed in Appendix lAl and [Bl 
respectively. In Sec. IIV1 we present a bath correlation 
function parametrization scheme that will be used for 
the later development of differential EOM formalism for 
general non-Markovian dissipation systems. In Sec.fVl 
we exploit a model dissipation system to illustrate the 
key ingredients in the IGF-COPI construction of the hi- 
erarchical EOM, augmented with residue correction and 
truncation. The final exact EOM formalism for general 
cases is presented in Sec. lVIl It is in principle appli- 
cable for arbitrary non-Markovian dissipation systems, 
interacting with the Fermionic/Bosonic grand canonical 
bath ensemble at any temperature, including T = 0. Fi- 



2 



nally, Sec. IVIII presents the comments, discussions, and 
concluding remarks of this work. 



II. TOTAL COMPOSITE HAMILTONIAN AND 
FLUCTUATION-DISSIPATION THEORY 

A. The description of system— bath coupling 

Let us start with the total system-plus-bath composite 
Hamiltonian, which in the stochastic bath /i B -interaction 
picture is given by 

H T (t) = H(t) + J2 [Wafl(t) + UtWl] . (1) 
a 

The deterministic Hamiltonian H(t) for the reduced sys- 
tem may also involve a time-dependent external coherent 
field drive. The second term in the right-hand-side (rhs) 
of Eq. ([1]) denotes the stochastic system-bath interaction. 
It is expressed in the multiple dissipative mode decom- 
position form, where {W a } and {f a (t) = e lflBt f a e- l,lBt } 
are the system and bath operators, respectively. They 
are in general non-Hermitian. 

The interaction bath operators {f a (t)} are also as- 
sumed to the Gaussian stochastic processes (a la the cen- 
tral limit theorem in statistics). The required Gaussian 
statistics is satisfied strictly when the individual f a (/J) is 
a linear combination of the annihilation (creation) oper- 
ators of the uncorrelated Fermion/Boson bath particles. 
The nonperturbative QDT formalisms to be developed 
later in this work will involve no further approximation. 

Throughout this work, we set h = 1 and the inverse 
temperature (3 = l/(k B T). Denote (0) B = tr B (Op B q ), 
with the grand canonical equilibrium density operator of 
the bare bath being given by 

-/3(/i B -AiiV) 
P^ = — — ■ (2) 

tr B [e-W'B-MJV)] 

Here, \i denotes the chemical potential of the bath. The 
particle number operator N — J2k c ! c fc °f the bath 
reservoirs commutes with the bare bath Hamiltonian, 
[N, h B ] — 0. In the case of a canonical bath ensemble, 
the number of particles is conserved and no longer a dy- 
namical variable. The resulting canonical p B q cx is 
equivalent to set the bath chemical potential fi = 0. 

For the later development, we also recast the system- 
bath coupling, i.e. the second term in the rhs of Eq. fTJ), 
as 

#'(*)= (3) 

a,cr 

Here, W% = (W^Y that can be either W a or W]; and 
similarly for f a = (/ CT )^. To be more specific, the bath 
operators {f a } considered in this work are linear com- 
binations of the annihilation operators, f a = ^2 k t a kCk, 



with CjC k ic^Cj = 5jk for the Fermion (+) or Boson (— ) 
particles in the grand canonical bath reservoirs. 

Note that as the stochastic bath operators f^(t) = 
[fa (*)] are the linear combinations of either the creation 
or the annihilation operators of the particles, we have 
(f%(t)} B = (fa{t)fb( T )) B = 0. These stochastic bath in- 
teraction operators are also assumed to follow the Wick's 
theorem for the thermodynamic Gaussian average that 
depends only on the two-time correlation functions, 

C^\t-r) = {fl{t)f b {T)) B , (4a) 
C { ab \t-T) = (f a (t)ft(r)) B . (4b) 

B. Fermion versus Boson bath reservoirs: 
Fluctuation-dissipation theorem 

Now we present the fluctuation-dissipation theorem 
(FDT), and the related symmetry and detailed-balance 
relations for the involving interaction bath correlation 
functions. The key steps to the final results will be de- 
tailed in Appendix [3] To avoid the confusion of later 
using ± to represent Fermion/Boson bath ensemble, let 
us recast Eqs. (0} as 

C i f b ) (t) = (fZ(t)ff(0)) B . (5) 

They satisfy the symmetry and the detailed-balance re- 
lations as [cf. Eqs. (fM|)-(|X3)] 

[CffW]* = Cfji-t) = elegit - t/3). (6) 
Their frequency-domain counterparts are 

C^H = [Cg(u>)]* = e^+^C^i-u), (7) 
with the corresponding spectral functions of 

*) = 2/ M <fc W *(^ (8) 
The spectrum positivity that reads 

C7&>M>0, Ctt{u)C%\u)>\C%{u)\\ (9) 

can be readily verified via the same method used in the 
Appendix A of Ref. Hit The symmetry and detailed- 
balance relations [Eqs. (|6|) or Eqs. (0] are general for a 
grand ensemble, no matter it is of Fermion or Boson in 
nature. This nature enters in terms of FDT via the spec- 
tral density functions as follows. 

Consider h B = J2k e k c l c k and f a = J^k^kCk- In this 
case, the interaction spectral density functions are de- 
fined physically for to > min(efc) = as 

J ab (u} > 0) = TT^2t ak t* bk S(uj - e k ) = J ba (u>). (10) 

k 

Mathematically, one can extend their definition to the 
lo < region, by setting J a b{^>) = ±Jba{—u) for the 
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Fermion/Boson (+/— ) particles for the reasons that will 
become self-evident soon. As a result, 

Jab(u) = Tr^2[t ak tl k 6{uj - e k ) ± t* ak t bk S(uj + e fc )]. (11) 

k 

On the other hand, 

[/aW,/ b t (0)] ± =E« fe e-^ t , (12a) 

k 

[P a (t)J b (0)} ± = ±J2t: k t bk e^\ (12b) 

fe 

assuming c-numbers in this case. Consequently, Eq. (jllj) 
can be recast as 

J ab (w) = - / ^ t ([/ a (t),/ b t (0)] ± ) B 

/"DO 

dte^([p a (t),f b (0)] ± ) B , (13) 

O 

or [cf. Eqs. Q and ©] 

^H = ci; ) M±(t ) H. (w) 

Note the symmetry relations here, 

■4a M = ±-4&(-<*>) = Jafe(^). (15) 
From the second identity of Eq. ([7]) , we have also 



r {a), s, -4bM 



or equivalently 



(16) 



(17) 



This is the FDT for the Fermion (+) or Boson (— ) grand 
canonical bath ensembles. 

In concluding this section, let us make some comments 
on Eq. (fl4|) . which can in fact be considered as the work- 
ing definition of the spectral density functions in general. 
For the linear bath interaction model, Eq. (|14p is equiva- 
lent to the conventional expression, Eq. (jTOJ) or Eq. (fTTj) . 
in which J ab (u>) is independent of temperature and chem- 
ical potential. As the working definition, Eq. (jT5J) may 
also be useful even for the cases where the bath interac- 
tion depends nonlinearly on the annihilation (creation) 
operators of the Fermion/Boson bath ensembles. In such 
cases, Eq. (|14l) that does not depend on the cr-index may 
involve with a certain mean-field approximation and the 
resulting J ab (uj) remains explicitly independent of tem- 
perature and chemical potential. The implicit depen- 
dence of Jab(w) on temperature and/or chemical poten- 
tial may however arise from the mean-field values of t ak 
or t bk , as in Eq. ([TO]) or (fTTj) . 

Consider now the value of J ab (uj = 0). As in- 
ferred from the second identity of Eq. (|15p. J a b(0) = 



±J* fc (0), implying that J a {,(0) is pure real/imaginary for 
a Fermion/Boson bath interaction. For the spectral den- 
sity function, assuming it is analytical at w = 0, the 
Boson bath interaction will also assume 



0; (Boson) 



j Qb (o) = [(i- e -^; ) H]_ 0;M=0 

as inferred from Eq. (fTB"]) and from the fact that J a h(w) 
does not depend on p explicitly. 

We are interested in nonperturbative QDT, in the 
presence of non-Markovian multi-mode interaction with 
grand canonical Fermion/Boson bath. Both the path 
integral and the differential EOM formalisms of QDT 
will be considered. The former will be presented in the 
coming section, while the latter will be constructed in 
Sec lVIi on the basis of calculus (time-derivative) on the 
exact path integral expressions. In order to have the 
EOM hierarchically coupled and closed for a general non- 
Markovian case, the bath correlation functions CjZ (t) 
should be parameterized in a proper form. The afore- 
mentioned properties of J b(w = 0) will be used there, to- 
gether with Eqs. (15]) and Eq. ([17]), to the FDT-preserved 

parametrization of C^ b (t) in Sec. IIVI that participates in 
the construction of the hierarchical EOM for general non- 
Markovian dissipation dynamics. 



III. PATH INTEGRAL FORMALISM AND 
DISSIPATION FUNCTIONALS 

In this section, we construct the path integral formal- 
ism of quantum dissipation. For the sake of clarity, we 
present only the final results here. The derivation of the 
path integral formalism will be given in Appendix [Bl It 
is similar to our previous work on the Bosonic canoni- 
cal bath interaction cases, where the system and bath 
operators, {W a } and {/ a }, were also assumed to be 
Hermitian. 5 -^ Briefly speaking, it employs the Wick's 
theorem for the thermodynamic Gaussian average, fol- 
lowed by using the symmetry and detailed-balance re- 
lations of Eq. ©. It involves also the ansatz of the 
initially factorizable total density operator, pt(^o) — 
p(to)p < B 1 , with p| q being given by Eq. ([2]). Note that 
when the initial time to — > — oo this ansatz imposes no 
approximation^^ 

Unlike its EOM counterpart that can be expressed 
in the operator level, the path integral formalism goes 
with a representation. Let {|a)} be a basis set in the 
system subspace, and a = (a, a') for abbreviation, so 
that p(ot,t) = p(a,a',t). Denote U(t,to) as the reduced 
Liouville-space propagator, p(t) = U(t,to)p(ta), which in 
the a-representation reads 

p(a,t) = J da U(a,t;a ,to)p{ao,to), (18) 

with the path integral expression of 11 
/■«[*] 

U(a,t;a ,t )= Va e lS[a] F[a} e - lS[a ] . (19) 

J ot a [ta\ 
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S [a] is the classical action functional of the reduced sys- 
tem, evaluated along a path a(r), with the constraints 
that two ending points a(io) = ao and a{t) = a are 
fixed. The key quantity in the path integral formalism is 
the Feynman- Vernon influence functional T[a\ that will 
be presented soon. 

In connection to the later development of EOM for- 
malism, we denote hereafter a = (aa') for a pair of dissi- 
pation modes and introduce 

f'C; {«}) = W£?(t; {«}) - W%\t; {a'}), (20) 



where [noting that W&' = W^J and Ca 



C<f}] 



W£\t; {a}) = / drC^(t - r)W^[a{r% (21a) 



in 



W^(t; {a'}) = / dr[C^(t - r)]*W^ [a'(r)].(21b) 
Denote also 

W Q >] = WMt)] - W^[a'(t)] = W a a - Wf. (22) 

It is in fact the -commutator in the path-integral 
representation, as it depends only at the fixed ending 
points of the path. The time-dependence here can thus 
be removed. 

The final expressions for the influence functional are 
[cf. Eq. l(B7]) ] 

T[a] = expj-J^ dr7e[T;{a}]| , (23a) 



with 



ft[i;{a}]=£w£[a]wW(t;{a}). 



(23b) 



Here, T is expressed in terms of its exponent time- 
integrand 1Z. It is in fact the time-local dissipation super- 
operator in the path integral representation, as the time 
derivative of Eq. (|23a[) . 



dtT 



-TIT, 



(24) 



leads to d t U = -iCU - K{t)U [cf. Eq. ([15])]. or equiva- 
lently p = -iC p - R(t)p [cf. Eq. Ip]l]. Therefore, 71 of 
Eq. (|23b[) can be termed as the dissipation functional. 

It is noticed that the chemical potential p does not ap- 
pear explicitly in the dissipation functional 1Z. The above 
exact path integral formalism, Eqs. (fl"8|) - ([23"]) . is formally 
the same for both the canonical and the grand canoni- 
cal bath interactions, no matter the bath is of Fermion 
or Boson, as long as the Gaussian bath statistics is ap- 
plicable. Apparently, in the case of W a = W\ = Q a 
and fa = fa = Fa, the sign-index a ( " + " and " — " ) 
can be omitted, and Eqs. ([20)) - (j2"3"]) recovers formally the 
Eqs. (10)-(13) of Ref. [5l|, where the Bosonic canonical 
bath ensembles were also assumed. The different nature 
of the quantum bath ensemble is embedded implicitly in 
the correlation functions C^\t), via the FDT and sym- 
metry relations as described in Sec. Ill Bl see Eqs. (fT5)) 
and (fT7|) . 



IV. NON-MARKOVIAN BATH VIA 
PARAMETERIZATION 

To construct a hierarchically coupled set of EOM for 
a general non-Markovian dissipation, C^a\t) should be 
expressed in a proper form, such as an exponential se- 
ries expansion ) 49 i 50 i 51 i 56 which shall also satisfy the FDT 
[Eq. (fI7|)]. To that end, let us consider the following 
parametrization form on the bath spectral density func- 
tions (J = J aa >), 



J a (u) = J5(u) + 



K 

E 

*:=i 



± 



( W - LOtf + ( 7fe °) 2 



(^ + ^)2 + ( 7 «)2 



with the Drude term 



^ + ( 7 g) 2 



or 



^ 2 + hS) 2 ' 



(25a) 



(25b) 



for the Fermion (+) or Boson (— ) bath ensembles, re- 
spectively. 

All involving parameters are real; o>™ and 7° (includ- 
ing 7™) are positive as well. They satisfy [cf. Eq. |[T5])] 

{uf Af Xf) = (Wfc' ,7fc'°, - Cfc' .-Cfc'°). The cor- 
responding bath correlation functions can then be ob- 
tained via the FDT of Eq. (fT7|) . using the contour inte- 
gration method. The final results read [M — > 00) 



2K 

C£\t) =vi a) e-^ t + Y,4 a Uj(t)+ £ ^n\t). (26) 

j — l m— 1 



A I 



» 



Va,j and Vj 



77^°j and similarly for 



Hereafter, r\ J 

other parameters that depend both (a) and a. 
The first two terms in the rhs of Eq. ([26]) . where 



0J fc (f)=cos(wJt)exp(->yj?t), 



(27a) 
(27b) 



arise from the poles of J a {z) in Eq. (I25|) . The involving 
^-coefficients are 



gj/5(7S+ cri i 1 ) -f 1 



3 i/3(7g+tri/j) _ ^ ' 



for the Fermion (+) or Boson (— ), respectively; JKfc— 1 = 
-tAgi?£ o) ± iA£i^ a) and = \%D^ ] ± Agz^ o) . Here, 

^° = m + csh + cm/pti?) = ckt, while 4 a) = 

{l±exp[(3(uj% + i~f% — ap)}}~ 1 , and J> i [ a ' 1 is similar as 
but with — instead. 

The last term of Eq. (j26|) arises from the poles of Mat- 
subara frequencies, 



7 m = (2m — l)ir/{3 or 2ra7r//3 



(28) 
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for the Fermion or Boson bath ensemble, respectively. 
CW(t) = fi^ exp airit] = T[C^(t)}*, (29) 



like Eq. (J3TJ) as 



n 



E 



(32) 



with 



1 



(a) = 



-i(2//3)J a (-ij m -a^ = T[&" 



(30) 



The last identity in each of the above two equations is 
originated from J a (—z) = ±[J (z)]*, the analytical con- 
tinuation of J a (-ui) = ±J*(uj) [cf. Eqs.fUlJ)]. It will 
show that this property leads to some simplification in 
treating the Matsubara-frequency contributions. In the 
canonical ensemble, the Matsubara coefficients fjm of 
Eq. (|30p . which do not depend on a and chemical poten- 
tial /j,, are real and pure imaginary for Boson and Fermion 
bath ensembles, respectively. 

To conclude this section, let us make some com- 
ments on the parameterized (t) of Eq. (j2l)|) , for which 
the hierarchical EOM formalism can be constructed via 
the IGF-COPI method 50 ' 51 without approximations (cf. 
Sec. lVI|) . In principle, Eq. (|26|) can be exact for a gen- 
eral non-Markovian dissipation, if the involving K and 
M are sufficiently large. The resulting EOM formalism 
is also exact; but its size grows in a power law as K or 
M increases. The exact evaluation of complex dissipa- 
tion would rapidly become extremely tedious. In prac- 
tise, the K and M in Eq. (|26[) have to be finite in the 
construction of the EOM formalism. To complete the 
theory, the residue correction due to the difference be- 
tween the true and the parameterized correlation func- 
tions, SC^ = ci°cxa — C^a \ should also be incorporated 
into the final EOM formalism. 



APPROACH TO THE EXACT EQUATIONS 
OF MOTION FORMALISM: PRINCIPLES 



A. Hierarchy construction 

For the sake of clarity, we exploit the broadband dis- 
sipation case to illustrate the IGF-COPI approach to 
the hierarchical EOM, together with the residue cor- 
rection and hierarchy truncation that will be demon- 
strated in the next subsection. In the broadband 
limit, exp(— j"t) — > t?d 7d<^(*) an d the 77^ — term 
in Eq. (|2"6")) can be completely neglected. We will show in 
Sec. lVBl that the <5(i)-like component can be easily taken 
into account in terms of the residue correction. 

Considered here explicitly for the illustration of hier- 
archy construction are only the Matsubara terms, 



M 



C { °\t) = C%\t) = J2 vtte-^-^. (31) 



m — 1 



m— 1 



The dissipative functional 7Z [Eq. (|23b[) ] , as it appears lin- 
early in the bath correlation functions, can be expanded 



Here, wffl are defined by Eqs. igD]) and flST) but with the 

bath correlation functions being replaced by Cm . We 
shall hereafter omit the explicit path-integral variables 
dependence whenever it does not cause confusion. 

Note that the Matsubara terms possess two special 
properties: C${t) = ±[c£\t)]* [cf. Eq. @] and the 
fact that their time constants (ff m — ai/i) are dissipation- 
mode independent. As a result, Eq. (|3"2")) can have the 
summation over a' performed, resulting in 



(33) 



The involving Eqs. (|2"D|) and (12"T]) are rearranged for the 
Matsubara terms as 



Ha) 



(w {a) 



(34) 



with 



W$ = fdr 6$ (t - r)W° a , [a(r)} , (35a) 

a' Jt 

Wl a) = E Ar<?W(t - r)C-K(r)]. (35b) 



We have 



where 



dt2g> = - (j m - ai[i)Z$ , (36) 



dg> = -iJ2n£HWZ±Wp)- (37) 



We are now in the position to construct the EOM 
via the calculus (time-derivative) on the path integral 
formalismi 50 i 51 The time derivative on the action func- 
tional parts contributes to the coherent dynamics of 
—iCU, and thus can be included into the final EOM 
formalism. However, the time derivative of the influ- 
ence functional, Eq. (|24[) . will generate other auxiliary 
influence junctionals (AIFs). This process goes on pro- 
gressively and hierarchically as the time derivatives on 
AIFs is continued. It is also noticed that time deriva- 
tive of individual Zm is closed by itself [Eq. ([55]) ]. We 
will see soon that this property, together with their ex- 
clusive relation to the dissipation functional as Eq. (|33p . 

make {zffl} the complete set of the IGFs for the desired 
hierarchy construction. 

Let us start with the time derivative on the primary 
influence functional, Eq. (|24|) . In the present demonstra- 
tive case, it reads 



(38) 
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Generated here are the first-tier AIFs, defined as 

J*$ = &$T. (39) 

Compared with the primary J- [Eq. (|23j) ] whose leading 

term is 1, the first-tier AIFs {Tm^} are of the second- 
order in the system-bath coupling as their leading con- 
tributions. 

To proceed, the time derivatives on those first-tier 
AIFs are also needed. We obtain [cf. Eqs. ([3"S|) - (j3"9"]) ], 

d t T^ = \C^T - (j m - <ni*)J*$\ + Z^{d t T). (40) 



Here. 



(oo ) 



(41) 



with (denoting 



•'mm 1 — m J ~ m' X^m n m' ) ' 



(42) 



being the next tier (second-tier) AIFs. Apparently, the 

leading contributions in are of the fourth-order in 

the system-bath coupling. 

Continue on applying Eqs. (|3"o ]) -([3"8" ]) for Eq. (g^J) 

dtJtS = [d t (Z^Z^?)]T+ (ztfZ%?)(d t F). (43) 

The first term in the rhs above will lead to the depen- 
dence on the tier-minus-one AIFs, arising from the Cm 
term of Eq. (|36|) . while the (dt .F)— term above leads to the 
tier-plus-one AIFs. 

Apparently, those distinct {Zm^} are IGFs, since all 
involving AIFs can be generically expressed as 

■Tmm'-m" ~ [ Z m Z m> ' ' ' Z m" K' I 44 J 

In order to eliminate possible duplications due to the 
indexes permutation/repetition that result in the same 
AIF, we arrange all distinct {Z^} in a specified sequence 
and rearrange Eq. (|44|) accordingly as 



^{ n m^y 



(45) 



Here n = {ft™' 1 } is the index-set of nonnegative integers 
for the number of occurrences on the distinct zffl in- 
dividuals in the specified sequence. Note that the total 
number of nonnegative integers in the index-set n is 2Mq, 
where M is the number of Matsubara terms considered, 
q the number of dissipative mode a, while the factor 2 
arises from a = + , — . 

To specify the hierarchical dependence of d t T„ , we de- 

{... *(») 



note also the index-set nf , = 



, n m ' ± 1, • • • } that 



firn^ ± 1. The time derivative of T n [Eq. (|45|) ] can be 
carried out, again by using Eqs. ([3^)) and (f3"8"|) , resulting 
in 

d t Tn = -(% + in»)K + £ nt ] ^^- a) 

a.(T,m 



Here. 



Tn = E ft m ^ = E^U + ^mllm, (47a) 



= -E™i a) = E^" 



-l-n£l]. (47b) 



The auxiliary reduced density operators can now be 
defined as (including also p = po) 

p n (t)=U n (t,t )p{t ), (48) 

with the auxiliary propagators of [cf. Eq. (|19p ] 

U n (cx,t;a ,to)= / Va.e lS ^T n [a]e~ lS ^ a l (49) 

J a 

The hierarchically coupled Eqs. (|4^|) can now be recast as 
p n = - [i(C + np) + 7n] Pn + E ft m4 B V fi - ) 

a,(Tjn 

-i E w >*u, m - (50) 

a,a,m 

Here, Wf and which were given in the path integral 
representation as Eqs. (f2"2"|) and (|3"T|) , respectively, are now 
the reduced Liouville-space operators; i.e., 

w° a 6 = [wz,d], c^6 = -iX^W.fyf (5i) 

a' 

B. Residue correction and the hierarchy truncation 

We have illustrated the IGF-COPI approach to the 
hierarchical EOM formalism. The same method will be 
used in Sec. IVII to establish the EOM formalism for a gen- 
eral non-Markovian dissipation characterized by Eq. (j2l)|) 
without approximation. The final hierarchical EOM [cf. 
Eq. I|74p] can be expressed as 



Pn 



[i(C+fip)+-f n +STl(t)}p n +p^ } +p i n - i +p i „ +} . (52) 



deviates from n only by changing the specified hm to 



The hierarchy-down p„~ } and hierarchy-up pi +} are sim- 
ilar as the last two terms in Eq. (j50")l . but with the Drudc 
and the (^"-components of the bath correlation functions 
of Eq. (f2"B]) included. The hierarchical-swap pj^ } arises 
completely from the ^"-components, thus have no cor- 
respondence in Eq. 
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Involved in Eq. (|52|) is also the residue correction 81Z(t) 
to the dissipation. It accounts for the practical differ- 
ence between the exact Ca^Jjb) and the parameter- 
ized ci <T ' ) (t) by Eq. where K and M are finite. In 
the following we shall present the principle of residue 
correction^ to construct the residue dissipation 81Z, fol- 
lowed by an efficient scheme of hierarchy truncation. 

Let us first verify that the residue dissipation correc- 
tion 6TZ = ft cxa - K, due to 8C ( a a) = c£% a - c£\ 
does enter at the global level, participating in every tier 
of the hierarchy as Eq. (|52|) . As inferred from Eq. (|23al) . 
the exact influence functional of primary interest reads 
in path integral representation as J rcxa = J 7 F rcsi , with 
dtT = —1ZJ- and dtJ- Icsl = —51ZF rcs i. The auxiliary in- 
fluence functionals defined, for example in Eq. (|45|) , for 
the construction of the hierarchical EOM should now be 
replaced by ^-"° xa = J- u J- vcsl . Its time derivative is then 
9 t ^ xa = (dtT^T^si-SIlT^. Together with the estab- 
lished dtT n such as Eq. (|46|) before residue correction, the 
above expression verifies the fact that STZ does enter into 
every tier of the final hierarchical EOM, as in Eq. (|52[) . 

We are now in the position to present some practically 
useful expressions of 81Z(t). Note that its exact expres- 
sion in path integral representation is [cf. Eq. (|23b|) with 
Eqs. O and ([21]) ] 

6K = J2w!{(x}6W^(t;{ a }), (53) 

a. a 

where 

5WP = SW^(t; {a}) - SW^\t; {a'}), (54) 

with 

= fdr SC^ (t - t)W:, [a(r)], (55a) 

SW'W = fdr (t t)] * W° a , \ol{r)\ . (55b) 

In the Markovian-residue limit, Eqs. (fS"5"|) reduce to 

WMK6CW(t)WZ[a(t)], (56) 
and similar for SWa^ ; with 




Note that the system variable WZ in Eqs. (fB"6")) is now 
evaluated at the ending time of the path integral where 
a(t) = a and ex! if) = a 1 are fixed. As results, Eq. (|53|) 
for the Markovian-residue dissipation can be expressed 
in the operator level as 

8KO = 5C£Xt)WZ,6- [5&iXt)]*dWZ\. (58) 



Alternatively, one can exploit various second-order 
QDT formalism^ for the residue 5R via, for example, 

SW<?\t) « fdrSC^it - T)e~ iC ^W:,. (59) 

The dissipation-free propagator (assuming also time- 
independent system Hamiltonian for clarity) is used here 
to connect W", [a(r)] in Eq. (|55|) to its value at the fixed 
ending path point of a(t) = a. It is the reason that the 
a-representation has been removed from Eq. (|59p , as it is 

now treated as an operator; the same for 5W'^ ] . There- 
fore, the weak-residue 81Z assumes the operator form as 

8KO = ^2[W!,5WtKt)6 -OSW^\t)}. (60) 

a.rr 

Apparently, the Markovian-residue limit is the special 
case of the weak- residue treatment. 

To the end of this section, let us show that the prin- 
ciple of residue correction presented above can also be 
used for the hierarchy truncation. Recall the hierarchical 
AIF considered in the previous subsection, T„ of Eq. (|45|) . 
Implied there is a trivial identity that 

Note that and W^ a) defined in Eqs. (J35D as their 

path integral expressions are of the same mathemati- 
cal structure as Eq. ([55]) . Therefore, the afore-described 
method, via either the Markovian or second-order ap- 
proximation scheme to the explicit operator form of 
can be adopted locally at the desired anchoring 
level N. Thus, the truncation scheme of 

P N + )im = ^°Vn = ^[W^,pn]±, (62) 

in which Wm assumes its approximate operator form is 
established. The choice of anchoring index-set N will be 
specified at the concluding part of the coming section, 
where the hierarchical EOM formalism for the general 
non-Markovian dissipation is treated. 

VI. EQUATIONS OF MOTION THEORY FOR 
GENERAL NON-MARKOVIAN DISSIPATION 

A. Influence generating functionals 

We are now in the position to the IGF-COPI con- 
struction of the exact hierarchical EOM formalism for 
the parameterized of Eq. (j2"6")l . The corresponding 
dissipation functional can be expressed as [cf. Eqs. (|23b[) 
and (021)] 

n = i £ w:4 a) + i E w « [4 a) + n (a) ] 

a, a a,a,k 

+i £ WfzW. (63) 

a,fT,m 



8 



The last term arising from the Matsubara contribution 
had been treated in detail in Sec. IV A[ see Eqs. (j3"3"|) - (f3"7) . 
Apparently, all these composite {X, Y, Z}-functionals are 
IGFs; but whether they are completed or not should be 
checked with their time derivatives; see Sec. IV Al 

Consider the Drude-IGFs in Eq. (|63|) , they are given 

by 



Z< a) = -i 



i(vt } M a) -V^*W^), (64a) 



These six classes of (XKZ)-functionals constitute now a 
complete set of IGFs for the hierarchy construction. The 
general expression for all auxiliary influence functionals 
in the hierarchy can then be expressed as 

= { II [rt a) f Ll (4 a) f* Fi a) f Ll (4 a) f* 

a.k 



(72) 



with 



Jtn 



)}. (64b) 



We obtain 



(a) 



Here 



(65) 



(66) 



Therefore, {Z^} constitute the complete IGFs as the 
Drude-term contribution is concerned. 

Turn now to the {X, Y"}-IGFs. They are given by 



(67a) 



yt a) ^ V^W^\), (67b) 



with [see Eq. (gTJ) for tf>(t)] 



The index in T u is specified by the involving nonnegative 
integers, 

n = (n^\ 4f, n { S\ ; j = 1, • -,2K; m = l,- ,m) . (73) 

Therefore, the total number of the nonnegative integers 
in the index-set n is 2(p + 4Kp + Mq), where q denotes 
the number of dissipative modes, and p < q 2 the number 
of nonzero dissipative mode pairs. In terms of Eqs. (|48[) 
and (|49|) . the IGF-COPI approach to the hierarchical 
EOM is now completed by considering the time derivative 
of T u [Eq. (l72|) ]. which can be readily carried out using 
Eqs. HMD, flU, dnSD and (JTUD ■ 



B. Hierarchical equations of motion 

The final hierarchical EOM for the associating auxil- 
iary density operators read 



p n = ~[i(C + hfi) + 7n]/5n + pn + Pn ' + Pn 



•>{+> 



(74) 



w. 



(a) 



dT^(t-r)W:,[a(r)}. 



(68) 



in 



Unlike the Z-functionals for the Drude and Matsub- 
ara components, the time derivatives of the X- and Y- 
functionals are closed, together with two additional func- 
tionals, given by 



i^W^U-^W'^), (69a) 



Y, 



(a) _ 



iiv^W^-v^W^). (69b) 



We have 



and 



d t X { k a) = A[ a) - tfX<°> - ^ k X { k a) , (70a) 
d t X { k a) =^xi a) -^xi a) ; (70b) 



The /i-term in Eq. (|74[) arises from the corresponding con- 
tribution of Matsubara term [Eq. ([36])] given by Eq. (147b|) . 
The 7-terms in Eq. (174|) arise from the damping-terms of 
Eqs. (|70|). The resulting damping constant is given by 



n 2k * ,l 2k-l 



-So) 
l 2k~l 



Ik 



a,k 



(75) 



The second term in Eq. (|74|) stems from the (off- 
diagonal) swap-terms of Eqs. ([70)1 . It reads 



= J2 W fc("2fc-lPn- 



(o) 
U 2k Pn- , 



- n S*-lP"{; ) ,»_ 1 + "2feVn- 9t ) ■ (76) 



Here 



8 t Y k ia] = B[ a) - ^Y k (a) - w%Y k {a) , (70c) 



dtY k 



(a) 



(a) 



(a) 



4 Q) - - 

B k a) = -i 



(70d) 



(^WZ-^W'j), (71a) 
(rftLM-^WZ). (71b) 



The index-set n7^ . and nj~* . differ from n of Eq. (173 



only at the specified {(a), j} by 



'(a) J 



(n 



(°) ^( a )\ 



(n 



(a) 



1, n 



(a) 



1), 



(n( o) +l,fi< >-l)^n<°U$ o) ). 



The third term in Eq. (O stems from the (A, B, C)- 
terms of Eqs. (|70p . while the last term is from Eq. 
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They are the hierarchy-down and hierarchy-up contribu- 
tions, respectively, and given by 

a.a,k 

and (i=D,l,--- ,2K) 

The reduced Liouville-space operator in Eq. (175)) 
and Cm ' involved in Eq. (fTT)) were given by Eq. (fSTj) . 

In Eq. JZB, A[ a \ B L k a \ and C^> a) denote the reduced 
Liouville-space operator counterparts of Eqs. (|66|) and 

CD 

4 Q) - -<(^W50 - vit OWl), (79a) 

S^O = ^{r^_ x W:,6 - r,^ W;,), (79b) 

d a) = -i (r^W^O - vi a) *OW^ , (79c) 

The index-set nf < . (n7 , . or n^, ) differs from n only 
by changing the specified or rim ] ) to nj a) ± 1 

(nj°'-l or nm'±l). Note that the (n^ l a ^) th -auxiliary re- 
duced density operators are not generated from Eq. (J78J) , 
since and Yfc do not appear in the influence phase in- 
tegrand 7Z [Eq. Q24|l]: they are rather generated from the 
EOM for the (n^ J ) th -auxiliary reduced density opera- 
tors via the involved swap {^}-terms [cf. Eq. (|76[)]. 

Note that, in Eq. (|74|) . the reduced density matrix 
Pq = p is of primary interest, and the initial conditions 
are p n {to) — p(ta)S n o, as inferred from their definitions. 
Setting the initial time to — > — oo, the established Hier- 
archical EOM is imposed no approximation. The initial 
conditions are now p(to) — 0, where to can be at any time 
before applying the the external time-dependent fields. 
The pulse-field induced dynamics will then be evaluated 
via Eq. (|74| in the presence of external field. 

C. Comments 

1. Residue dissipation 

The hierarchical EOM, Eqs. (J74])-(|79|), are exact for a 
general dissipation system that involves the parameter- 
ized bath correlation functions of Eqs. (j2l)|) . Due to the 
finite K and M of Eqs. ([26]) in practice, the complete the- 
ory can be established by considering the residue dissi- 
pation 51Z(t), due to the small difference between the ex- 
act and the parameterized (i) as described in section 
IVBl The final hierarchical EOM thus read [cf. Eq. ([52])] 

Pn = -[i{£+hp)+ ln +5K{t)} Pn +p^ } +p { n - } +p { n +i . (80) 



Actually, this hierarchical EOM Eq. (|80|) is valid for ar- 
bitrary temperature, including T = 0, at which the Mat- 
subara damping constant 7 m is meaningless and Mat- 
subara expansion is not valid. Based on the principle of 
the residue correction, considering the finite temperature 
in the Matsubara term, the difference between the exact 
(T = 0) and the parameterized ct } (t) (finite T) can be 
incorporated in the residue dissipation 8TZ(t). 

In the high temperature limit, a large Matsubara 
damping constant 7,„ oc kT leads to the Markovian-limit 
bath correlation of (t). In this case, the Matsubara- 
contribution can also be included in the residue dissipa- 
tion 51Z(t) as given by Eq. (|58|) which enter the theory 
globally at every hierarchical tier. 



2. Truncation consideration 

The hierarchy truncation via the principle of residue 
correction has also been demonstrated in subsection lVB] 
where only the Matsubara term is considered. In addition 
to the truncation of the Matsubara term of Eqs. (|6Tj) and 
(|62| , the same procedure is applied to the Drude and the 
4> terms (p N + and p N + ). 

'\ Q ),D 

Now we consider the choice of the anchor index-set N . 
The anchoring indexes can be specified for the individual 
constituents of the interaction bath correlation functions 
C [ a\t) of Eq. flU as 

NS, Ng, ^C; with k = 1, ■ ■ ■ , K + 1; m = 1, • • • , M. 

The closed set of hierarchically coupled EOM will then 
contain those p n , whose individual index-set consists of 
the nonnegative integers that are confined within 

e ct + + < a k + 4ti < N k, (81^) 

£.4> a) <A^ E CT ^<^ (81b) 

The anchor index-set N in p^ can now be defined as those 
with at least one upper limit being reached. The trun- 
cation can thus be made based on the formally exact 
relations such as Eqs. (|6Tj) and ([62]). 



VII. CONCLUDING REMARKS 

In summary, we have generalized the exact non- 
Markovian quantum dissipation theory , 50 ! 51 on the basis 
of the calculus-on-path-integral algorithm, to non- hermit 
coupling and non-Markovian grand canonical bath en- 
sembles (both for Fermion and Boson statistics). The 
resulting dissipation functionals path-integral expression 
is found to be of the same as the canonical bath with 
the assumption of they all satisfying the Gaussian statis- 
tics. The distinct fluctuation dissipation properties 
lead to different hierarchical EOM coupling structures. 
However, the difference appears only in the Matsubara 
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contributions, Boson-commutator versus the Fermion- 
anticommutator relations, due to the particularity of the 
Matsubara coefficients [Eq. ([30]) ]. In the high tempera- 
ture limit, the aforementioned distinction diminishes, as 
both the Fermion and Boson bath ensembles approach to 
the Boltzmann statistics. 

The hierarchical EOM formalism, which is in princi- 
ple equivalent to its path-integral counterpart, may be 
relatively tractable; see the Summary of Ref. [5l|. Its nu- 
merical implementation remains a challenge, especially 
for complicated, strong dissipation systems at the ex- 
tremely low temperature regime. Note that the hierarchy 
is not just for system-bath coupling strength, but more 
importantly also for bath correlation time scaled The 
proper choice of system to effectively reduce the non- 
Markovian system-bath coupling strength in Eq. is 
therefore still important both physically and numerically 
here. The standard approaches in this regard include the 
variation-principle based canonical transformation (e.g. 
the polaron picture) and the primary-bath-mode method 
(e.g. the solvation coordinate picture) i 57 ' 58 i 59 ' 60 i 61 

The present theory can account for the particle trans- 
fer between the system and the bath. Thus the prob- 
lems of quantum transport through molecular wires and 
electron transfer in complex molecular systems can be 
demonstrated by this theory. For instance, if W a (WJ) 
represents annihilation (creation) operator of the system 
in the coupling Hamiltonian [Eq. ([5])]. the measurable 
quantity, such as the current can be exactly expressed 
as (/(<)> - -i([N,H T (t)]) r = -iTr a J2 a {W aP Ut) - 
W^p a (t)}. Here p a (t) [Pa{t)] is the summation over the 
reduced density matrixes in the first-tier described by 
Eq. (|74[) with n = 1. This detail work will be treated 
elsewhere. The present theory can be further generalized 
to the system coupled with many reservoirs (including 
phonon bath). The work along this line is in progress 
and will be published elsewhere. 
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APPENDIX A: QUANTUM STATISTICAL 
PROPERTIES OF THE INTERACTION BATH 
CORRELATION FUNCTIONS: DERIVATION 

The symmetry relation - the first identity of Eq. ([6]). 
This can be obtained immediately from the definition 
of correlation functions, Eq. ([5]), together with the time- 
translation invariance and the trace cyclic invariance. 



The details are as follows [noting that /£ = (/£)t]. 

[C$>(t)]* ee tr B [e^e-^7 6 V B T 
= tr B [p c B V- e ^*/ a ff e- l/lBt ] 

= tv B [f?(o)f!(tW] 

=ci:\-t). (ai) 

The detailed-balance relation - the second identity of 
Eq. ©. Its derivation exploits the property of 

e*"/;r^ = e^tf. (A2) 

This relation arises from the linearity of f a = J^. t a jCj, 

together with the notions that /~ = f a and /+ = /J. 

Let us start with the derivation of Eq. (|A2|) . It is car- 
ried out by recognizing some basic quantum mechanics 
relations as follows. The first one is 

oo 

e A Be- A = e A B = —A n B, (A3) 

n=0 

where AB = [A, B]. We have therefore 

ftffo-ert = NT ^f>v;, (A4) 

whereA7 a CT = [A>,/ Q CT ]. 

We shall also use the quantum mechanics relations of 

[AB,C]=A[B,C] ± t[A,C] ± B, (A5) 

where [•,•]- = [v] an d [v]+ = {v} denote the com- 
mutator and anticommutator, respectively. By noticing 
that [cfc,cj]± = Skj and [cfc,Cj]± = for Fermions (+) or 

Bosons (— ), we have [N,cj] = Y^ki c k c k> c j] = — c j an d 
[N,c]] = c). Thus, [Nj a ] = -f a and [N, /]] = /t; i.e., 

AffZ = [N, f%] = aj a a . It immediately leads Eq. (|A4| to 
Eq. (|A2|) . which also implies 

efl(fcB-M*)/|( t ) e -/J(fcB-M*) = e al3tl f°(t - i/3), (A6) 

since [h B ,N] = and e /3/lB /^(t) e - /3,lB 

The derivation of the detailed-balance relation, the sec- 
ond identity of Eq. ([6]), is now straightforward, by using 
the last identity of Eq. (|A1[) , the trace-cyclic invariance, 
Eq. m and Eq. ££!]). We have 

cL CT) (-0=tr B [/f(tKV^(0)] 

= ti B [f!(t)e-^-^\f^(0)pl^ h --^} 

= ^ e 0(h B -»N)j-9^ e -P(h B -»N)f*( Q ^ 

= e^C s ab {t - »/3). (A7) 

The derivation of Eqs. ([6]), the symmetry and the 
detailed-balance relations in the time domain, is 
now completed. Their frequency-domain counterparts, 
Eqs. ([7]) with Eq. ([5]) , are then followed immediately. 
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APPENDIX B: PATH INTEGRAL FORMALISM: 
DERIVATION 

Let £/t(£, to; {f% (t)}) be the propagator of the total 
system and bath space, satisfying the Schrodinger equa- 
tion of the total Hamiltonian of Eq. ([1]) , 

dt U T = -i[H(t) + W s a r a {t)} U T . (Bl) 

a,cr 

The total density operator at time t is given by 

PT (t) = U T (t,t ;{fZ(t)})p T (t )U{(t,t ;{fZ,(t)}). (B2) 
Assuming it initially at time to is [cf. Eq. @] 

Pt(*o) = p(to) P c B q oc p{to)e-^-^\ (B3) 
the reduced density matrix of primary interest is then 
p(t)=ti B [p T (t)} 

= tx B [U T (t,t ;{f:(t)})p(t )e-^-^ 

t4(«,*b;{/5(«)})^e^»-"*)] 

'u T (t,to;{e^fZ(t-i[3)})p(to) 

= U(t,t )p(t ). (B4) 

The third identity is obtained by using the trace cyclic 
invariance and the relations of Eq. (|A6[) . 

Consider now the path integral expression of the re- 
duced Liouville-space propagator U(t,t ), as defined in 
the last identity of Eq. (|B4j) . Let {\a)} be a basis set in 
the system subspace. We have 

U T (a,t; a , to; {/*(*)}) 

= /pa^WexpJ-^ fw:[ a (r)]/;(r)}, 

V\{ol ',t;a ,t ;{fS>(t)}) 

C Vol e~ iS M exp U £ r dT ^ [ a '( r) ] f a , ( r) } . 

The action functional S[a] is related to the reduced sys- 
tem Hamiltonian H(t) only. The system operators {W 7 ^} 
have also been represented in path integral representa- 
tion. On the other hand, the stochastic bath dynamics 
variables {/^} remain in the operator form that involve 
in the time-ordered exponential functions. 

Substituting the above path integral expressions for 
Eq. (|B4|) leads to Eq. (fill)) , with the influence functional 
there being given by 

T = (cxp + {-»]T f drW°[a(T)]e^fZ(T - 0)} 



a,<r- t( > 

ft 



x exp_ 



drW!,[a'(rm(r)\) 

a',<j Jt ° 



For {/„ (t)} satisfying Gaussian statistics, the bath en- 
semble average in Eq. (|B5I) can be evaluated exactly by 
using the second-order cumulant expansion method, as 
the higher order cumulants are all vanished. This prop- 
erty, together with Eq. ([6]), leads to the expression 



$[a, a 1 ] 

= E fdT2rdr 1 {w a [a{T 2 )]Wl\a{T 1 )]C^J{T 2 - n ) 

a,a> Jt J to 

+Wl[a{T 2 )]W a ,[a{T 1 )]C { J{T 2 - n) 
+Wl,[a'{T 1 )]W a [a'{T 2 )]C { aa }*{r 2 -r 1 ) 
+W al [a'(r 1 )]W^(r 2 )}C { a + J*(r 2 - n)} 

E A 7 * fdrx{w a [a(r 2 )}Wl [a'(n)]C' >2 - n ) 
+Wl[a{T 2 )]W a ,W{r 1 )]C { +}*{T 2 n )] 



(ItK[t; {a}]. 



(B6) 



to 



for the influence functional exponent. 

In detail, we have used the symmetry relation [the 
first identity of Eq. |6])] in the third and fourth terms, 
and the detailed-balance relation [the second identity of 
Eq. |[S])] in the last two terms of the above expression. 
Some elementary algebra will lead to Eq. (|B6[) , in terms 
of 1Z = <9 t $, the expression, 



K[t;{a}}=J2({W a [a(t)]-W a [a'{t)]} 



x{W^[t;{a}]-W^ a V[t;{a'}]} 
+ {Wl[a{r)]-Wl[a'{r)]} 

x{W^)[f,M]-W^[t;{a'}}} 
£ {WMt)]- W s a [a\t)]} 



{W^[t;{a}]-Wy[t;{a'}}}, (B7) 



exp {— $[a, a']} 



(B5) 



with W$ = Wi° ] and W' a ^? = W^ a) being given by 
Eqs. dHJ). Using the definitions of Eqs. (gHD and ([22]). the 
above expressions can be recast as Eqs. 
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